

Polyfem Examples¶

This is a jupyter notebook. The "real" notebook can be found here.

Some imports for plotting

In [1]:

import plotly.offline as plotly
import plotly.graph_objs as go
import plotly.figure_factory as ff

#Necessary for the notebook
plotly.init_notebook_mode(connected=True)

algebra

In [2]:

import numpy as np

stuff for the animation

In [3]:

import ipywidgets as widgets
from ipywidgets import interact

and finallypolyfempy

In [4]:

import polyfempy as pf

Plotting utility¶

This code is not particularly interesting.

It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.

In [5]:

def create_plot_mesh_and_layout(vertices, connectivity, function):
    x = vertices[:,0]
    y = vertices[:,1]

    if vertices.shape[1] == 3:
        z = vertices[:,2]
    else:
        z = np.zeros(x.shape, dtype=x.dtype)

    if connectivity.shape[1] == 3:
        f = connectivity
    else:
        #Convert a tet-mesh into face based triangles.
        f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
        for i in range(len(connectivity)):
            f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
            f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
            f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
            f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])

    mesh = go.Mesh3d(x=x, y=y, z=z,
                     i=f[:,0], j=f[:,1], k=f[:,2],
                     intensity=function, flatshading=function is not None,
                     colorscale='Viridis',
                     contour=dict(show=True, color='#fff'),
                     hoverinfo="all")
    layout = go.Layout(scene=go.layout.Scene(
        aspectmode='data',

        xaxis=dict(
            autorange=True,
            showgrid=False,
            zeroline=False,
            showline=False,
            ticks='',
            showticklabels=False
        ),
        yaxis=dict(
            autorange=True,
            showgrid=False,
            zeroline=False,
            showline=False,
            ticks='',
            showticklabels=False
        ),
        zaxis=dict(
            autorange=True,
            showgrid=False,
            zeroline=False,
            showline=False,
            ticks='',
            showticklabels=False
        )
    ))

    return mesh, layout

In [6]:

def plot(vertices, connectivity, function, camera=None):
    mesh, layout = create_plot_mesh_and_layout(vertices, connectivity, function)
    if camera is not None:
        layout.scene.camera = camera

    fig = go.Figure(data=[mesh], layout=layout)

    plotly.iplot(fig)

Creates a quad mesh of n_pts x n_pts in the form of a regular grid

In [7]:

def create_quad_mesh(n_pts):
    extend = np.linspace(0,1,n_pts)
    x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy')
    pts = np.column_stack((x.ravel(), y.ravel()))

    faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32)

    index = 0
    for i in range(n_pts-1):
        for j in range(n_pts-1):
            faces[index, :] = np.array([
                j + i * n_pts,
                j+1 + i * n_pts,
                j+1 + (i+1) * n_pts,
                j + (i+1) * n_pts
            ])
            index = index + 1

    return pts, faces

Plate hole example¶¶

This is the python version of the plate with hole example explained here.

Set the mesh path

In [8]:

mesh_path = "plane_hole.obj"

create settings

In [9]:

settings = pf.Settings()

pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$

In [10]:

settings.discr_order = 1

normalize the mesh to be in $[0,1]^2$

In [11]:

settings.normalize_mesh = True

and choose Young's modulus and poisson ratio

In [12]:

settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)

We are use a linear material model

In [13]:

settings.tensor_formulation = pf.TensorFormulations.LinearElasticity

Next we setup the problem

In [14]:

problem = pf.GenericTensor()

sideset 1 has zero displacement in $x$

In [15]:

problem.set_displacement(1, [0, 0], [True, False])

sideset 4 has zero displacement in $y$

In [16]:

problem.set_displacement(4, [0, 0], [False, True])

sideset 3 has a force of [100, 0] applied

In [17]:

problem.set_force(3, [100, 0])

fianally set the problem

In [18]:

settings.set_problem(problem)

Solve!

In [70]:

solver = pf.Solver()

solver.settings(str(settings))
solver.load_mesh_from_path(mesh_path)

solver.solve()

Get the solution

In [20]:

[pts, tets, disp] = solver.get_sampled_solution()

diplace the mesh

In [21]:

vertices = pts + disp

and get the stresses

In [22]:

mises, _ = solver.get_sampled_mises_avg()

finally plot with the above code

In [23]:

top_camera = dict(
    up=dict(x=0, y=1, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=0, y=0, z=0.8)
)

plot(vertices, tets, mises, top_camera)

Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call

In [24]:

mises = solver.get_sampled_mises()
plot(vertices, tets, mises, top_camera)

Torsion¶

Non-linear example. We want to torque a 3D bar around the $z$ direction.

The example is really similar as the one just above.

Sets the mesh and create the settings. As before

In [25]:

mesh_path = "square_beam_h.HYBRID"

settings = pf.Settings()

It is an hex-mesh so we are using $Q_1$

In [26]:

settings.discr_order = 1

In this case the mesh has already the correct size.

In [27]:

settings.normalize_mesh = False

Choose Young's modulus and poisson ratio, as before

In [28]:

settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)

Differently from before we want a non-linear material model: NeoHookean

In [29]:

settings.tensor_formulation = pf.TensorFormulations.NeoHookean

and we want to do 5 steps of incremental loading to avoid ambiguities in the rotation

In [30]:

settings.nl_solver_rhs_steps = 5

Now we setup problem with fixed sideset, rotating an number of tours

In [31]:

problem = pf.Torsion()

sideset 5 is fixed

In [32]:

problem.fixed_boundary = 5

sideset 6 rotates

In [33]:

problem.turning_boundary = 6

around the $z$-axis (2)

In [34]:

problem.axis_coordiante = 2

by half a tour

In [35]:

problem.n_turns = 0.5

Now we choose a coarse visualization mesh

In [36]:

settings.vismesh_rel_area = 0.001

and set the problem and solve

In [37]:

settings.set_problem(problem)

#solving!
solver = pf.Solver()

solver.settings(str(settings))
solver.load_mesh_from_path(mesh_path)

solver.solve()

takes approx 1 min because it is a complicated non-linear problem!

Get solution and stesses as before

Since we want to show only the surface there is no need of getting the whole volume, so we set boundary_only to True

In [38]:

[pts, tets, disp] = solver.get_sampled_solution(boundary_only=True)
vertices = pts + disp
mises, _ = solver.get_sampled_mises_avg(boundary_only=True)

and plot the 3d result!

In [39]:

plot(vertices, tets, mises)

Fluid simulation¶

Create the mesh using the utility function

In [40]:

pts, faces = create_quad_mesh(50)

create settings

In [41]:

settings = pf.Settings()

settings.vismesh_rel_area = 0.001

pick linear $Q_2$ elements for velocity and $Q_1$ for pressure

In [42]:

settings.discr_order = 2
settings.pressure_discr_order = 1

Set the viscosity of the fluid

In [43]:

settings.set_material_params("viscosity", 1)

We select stokes as material model

In [44]:

settings.tensor_formulation = pf.TensorFormulations.Stokes

The default solver do not support mixed formulation, we need to choose UmfPackLU

In [45]:

settings.set_advanced_option("solver_type", "Eigen::UmfPackLU")

We use the standard Driven Cavity problem

In [46]:

problem = pf.DrivenCavity()

we set the problem

In [47]:

settings.set_problem(problem)

We create the solver and set the settings

In [48]:

solver = pf.Solver()
solver.settings(str(settings))

This time we are using pts and faces instead of loading from the disk

In [49]:

solver.set_mesh(pts, faces)

Solve!

In [50]:

solver.solve()

We now get the solution and the pressure

In [51]:

[pts, tris, velocity] = solver.get_sampled_solution()

and plot it!

In [52]:

top_camera = dict(
    up=dict(x=0, y=1, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=0, y=0, z=1.2)
)

plot(pts, tris, np.linalg.norm(velocity, axis=1), top_camera)

In [53]:

n_pts = len(pts)
scaling = 3
quiver = ff.create_quiver(
    pts[0:n_pts:10,0], pts[0:n_pts:10,1],
    scaling*velocity[0:n_pts:10,0], scaling*velocity[0:n_pts:10,1])

quiver.layout.yaxis.scaleanchor="x"
quiver.layout.yaxis.scaleratio=1
plotly.iplot(quiver)

Time dependent simulation¶

Create the mesh using the utility function

In [54]:

pts, faces = create_quad_mesh(50)

create settings

In [55]:

settings = pf.Settings()

pick linear $Q_1$ elements.

In [56]:

settings.discr_order = 1

and choose Young's modulus and poisson ratio

In [57]:

settings.set_material_params("E", 1)
settings.set_material_params("nu", 0.3)

We are use a linear material model

In [58]:

settings.tensor_formulation = pf.TensorFormulations.LinearElasticity

For efficienty in the browser we select a coarse vis mesh

In [59]:

settings.vismesh_rel_area = 0.001

We simulate from 0 to 10s and 50 steps.

In [60]:

settings.tend = 10
settings.time_steps = 50

Next we setup the problem, this doesnt have any parameters. It will...

In [61]:

problem = pf.Gravity()

we set the problem

In [62]:

settings.set_problem(problem)

We create the solver and set the settings

In [63]:

solver = pf.Solver()
solver.settings(str(settings))

This time we are using pts and faces instead of loading from the disk

In [64]:

solver.set_mesh(pts, faces)

Solve!

In [65]:

solver.solve()

Get the solution and the mises

In [66]:

pts = solver.get_sampled_points_frames()
tris = solver.get_sampled_connectivity_frames()
disp = solver.get_sampled_solution_frames()
mises = solver.get_sampled_mises_avg_frames()

Now we are ready to do the animation

First create a figure widget

In [67]:

top_camera = dict(
    up=dict(x=1, y=0, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=0, y=0, z=1.2)
)

mesh, layout = create_plot_mesh_and_layout(pts[-1],tris[-1],mises[-1])
mesh.cmin = 0
mesh.cmax = 0.25
layout = go.Layout(scene=go.layout.Scene(
        camera=top_camera,

        aspectratio = dict( x=1, y=1, z=1 ),
        aspectmode = 'manual',

        xaxis = dict(range=[-0.1, 1.1]),
        yaxis = dict(range=[-0.1, 1])
))
fig = go.FigureWidget(data=[mesh], layout=layout)

then we need the callback

In [68]:

def on_value_change(value):
    frame_index = value['new']

    fig.data[0].x=pts[frame_index][:,0] + disp[frame_index][:,0]
    fig.data[0].y=pts[frame_index][:,1] + disp[frame_index][:,1]

    fig.data[0].intensity = mises[frame_index]

finally we create an animation widged

In [69]:

play = widgets.Play(
    value=0,
    min=0,
    max=len(mises)-1,
    step=1,
    description="Press play",
    disabled=False
)

slider = widgets.IntSlider(min=0, max=len(mises)-1)

slider.observe(on_value_change, 'value')
play.observe(play, 'value')

widgets.jslink((play, 'value'), (slider, 'value'))
widgets.VBox((widgets.HBox((play, slider)), fig))

</div>